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ABSTRACT 

We investigate baroclinic instability in flow conditions relevant to hot extrasolar planets. The 
instability is important for transporting and mixing heat, as well as for influencing large-scale 
variability on the planets. Both linear normal mode analysis and non-linear initial value cal- 
culations are carried out - focusing on the freely-evolving, adiabatic situation. Using a high- 
resolution general circulation model (GCM) which solves the traditional primitive equations, 
we show that large-scale jets similar to those observed in current GCM simulations of hot ex- 
trasolar giant planets are likely to be baroclinically unstable on a timescale of few to few tens 
of planetary rotations, generating cyclones and anticyclones that drive weather systems. The 
growth rate and scale of the most unstable mode obtained in the linear analysis are in qual- 
itative, good agreement with the full non-linear calculations. In general, unstable jets evolve 
differently depending on their signs (eastward or westward), due to the change in sign of the 
jet curvature. For jets located at or near the equator, instability is strong at the flanks - but not 
at the core. Crucially, the instability is either poorly or not at all captured in simulations with 
low resolution and/or high artificial viscosity. Hence, the instability has not been observed or 
emphasized in past circulation studies of hot extrasolar planets. 

Key words: hydrodynamics - planets and satellites: atmospheres - methods: numerical - 
instabilities - turbulence - waves 



1 INTRODUCTION 

Baroclinic instability is a generic flow instability that occurs in ro- 
tating, stably- stratified fluids subject to a meridional (northward) 
temperature gradient. Examples of such a fluid are planetary atmo- 
spheres and oceans. The temperature gradient induces a vertical (al- 
titudinal) shear in the mean flow by thermal wind balance (e.g. Ped- 
losky 1987); hence, baroclinic flows are those that nominally vary 
in the vertical direction. The instability itself is important because 
it gives rise to large- and meso- scale weather systems on planets. It 
also serves as a source of turbulence, which has been invoked as the 
initial condition in some simulations of extrasolar planets to gener- 
ate plausible initial jet profiles (e.g. Cho et al. 2003, 2008). More 
importantly, the instability is a source of spatio-temporal variability 
which could be observed remotely. 

Baroclinic instability on extrasolar planets has not been stud- 
ied thus far. In this work we perform a simple linear analysis of 
a horizontally uniform jet in vertical shear. We also use a highly- 
accurate pseudospectral general circulation model (GCM) which 
solves the hydrostatic primitive equations to study the non-linear 
evolution of a non-uniform, gradient- wind balanced jet on an extra- 
solar planet. The primitive equations solved govern the large-scale 
dynamics of planetary atmospheres (e.g. Holton 1992; see also Cho 
et al. 2008 for some discussion relevant to the present work). Here 
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the main focus is on close-in gaseous planets, as they remain the 
best studied type of extrasolar planets thus far. However, much of 
the findings and discussion apply to hot extrasolar planets in gen- 
eral - regardless of whether a solid boundary is present or the ra- 
diatively stable layer extends deeply into the planet. 

For concreteness, we present calculations for a model planet 
with physical parameters appropriate for the extrasolar giant planet 
HD209458b (Table 1). We focus on the stability of broad, high- 
speed zonal jets - positive (eastward) at the equator and negative 
(westward) at high latitude - under adiabatic (i.e. heating and cool- 
ing rates balanced in the net) situation. By 'broad' we mean a width 
of ~ Ld, where Ld is the Rossby deformation length (section 2.1), 
and 'high-speed' means the speed is ~ 1000 m s _1 at the core of the 
jet. Such jets are commonly produced in diabatically-forced GCM 
simulations of close-in extrasolar giant planets (e.g. Showman et al. 
2008; Rauscher & Menou 2010; Thrastarson & Cho 2010). A study 
of adiabatic behaviour is needed because it provides the necessary 
baseline for comparing the instability under forced conditions and 
because, in many circumstances, the produced jets are not main- 
tained by the applied thermal forcing (but some flow-modified ver- 
sion, away from the specified radiative equilibrium). 

Our basic approach in this work is to carefully study baroclinic 
instability in sufficient generality, without complicating the funda- 
mental process with details which are still uncertain for extrasolar 
planets. The primary aim here is three-fold: 1) to ascertain the im- 
portance of baroclinic instability as a generic process operating on 
extrasolar planets; 2) to gain a better understanding of the outputs 
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from current extrasolar planet GCM simulations, made difficult by 
the complexity of solving the primitive equations accurately; and, 
3) to explore fundamental issues in baroclinic instability that have 
received less emphasis in traditional geophysical fluid dynamics 
studies, due to the markedly different parameter regime of many 
extrasolar planets compared to that of the Earth. 

The overall plan of the paper is as follows. Section 2 presents 
linear stability analysis. Linear growth rates and phase speeds are 
calculated for the traditional primitive equations on the '/3-plane', 
a differentially rotating plane tangent to the surface of the planet 
at a given latitude. In section 3 we present the non-linear evolu- 
tion of the instability, obtained from full numerical simulations. 
This section also presents the description of the numerical model 
and setup, as well as the non-convergence of under-resolved and/or 
over-dissipated simulations. The foundation for baroclinic life- 
cycle study is also laid in this section; a detailed discussion of the 
phenomenon is presented elsewhere, as are of forced evolution and 
transient growth. Recapitulation and discussion are given in sec- 
tion 4. 



2 LINEAR THEORY 

2.1 Charney-Stern-Pedlosky Criteria 

Necessary conditions for instability exist. These may be derived 
directly from global conservation of pseudoenergy and are given 
in Charney & Stern (1962) and Pedlosky (1964). Hence, we shall 
not derive the conditions here but simply list them for the reader's 
convenience. The conditions play an important role in this work, 
particularly in understanding the setup of the nonlinear initial value 
problem (section 3). 

Let q = a (cu + 2 ft) • V L (p be the potential vorticity, where 
cu is the relative vorticity, ft is the planetary vorticity, a is the spe- 
cific volume (= where p is the density) and V L is a gradient 
operator acting on a materially conserved field ip, which may be 
a function of temperature and pressure (e.g. potential temperature 
or entropy). Additionally, let x, y and z be the zonal (i), merid- 
ional (j) and vertical (k) directions, respectively. Given the zonal 
flow, U = U(y, z) i, and the basic state potential vorticity Q(y, z) 
such that q(x,y,z) = Q + q'(x,y,z), one of the following neces- 
sary criteria must be met for the onset of instability: 

(i) dQ j dy and dU / dz are opposite signs at the upper 
boundary 

(ii) dQ I dy and dU / dz are same signs at the lower 
boundary 

(iii) dU / dz is the same sign at the upper and lower 
boundaries - a condition that is distinct from 
condition (i) or (ii), if dQ/dy = 

(iv) dQ j dy changes sign somewhere in the interior. 

Note that Q — Q(U) and the prime denotes deviation from the 
basic state. 

For realistic flow profiles studied in section 3, the instability 
criterion is normally satisfied through criterion (iv). In addition, 
criteria (i) and (ii) are also satisfied in most cases. These conditions 
are useful for assessing stability of any basic flow configuration. 
However, they do not provide quantitative information, such as the 
growth rates of unstable modes and phase speeds of waves/eddies 
generated by the instability. 



Table 1. Numerical parameter values for HD209458b. Here g is surface 
gravity; R p is equatorial radius; fl is rotation rate; 1Z is gas constant; c p is 
specific heat at constant pressure; H is scale height; p T is reference pres- 
sure; T eq is equilibrium temperature; U is characteristic flow speed; and, 
TV is Brunt- Vaisala' frequency. 



Parameter 


Value 


Units 


9 


9.8 


m s -2 




10 8 


m 


n 


2.1xl0- 5 


s- 1 


n 


3.5xl0 3 


Jkg- 1 K" 


Cp 


1.23 xlO 4 


Jkg-iK- 


H 


5.8xl0 5 


m 


Pr 


10 5 


Pa 


L e q 


1500 


K 


u 


1000 


m s _1 


N 


2xl0- 3 


s- 1 



For the Earth, the stability analysis is typically based on the 
quasi-geostrophic (QG) theory, in which small Rossby number Ro 
and order unity Burger number Bu are assumed (Charney 1947; 
Eady 1949; Phillips 1951). Given the characteristic flow speed U, 
Coriolis parameter /, horizontal length scale L, and the Rossby de- 
formation length scale Ld, Ro and Bu are defined Ro = U / (fL) 
and Bu — (Ld /L) 2 , respectively. Here Ld = NH/ f is the Rossby 
deformation length, where N is the Brunt- Vaisala frequency; H is 
the characteristic vertical scale; and, f — 2fL sin 0, where ft = \ ft\ 
is the planetary rotation rate and 4> is the latitude. Note that both 
Ro and Bu vary with latitude. For example, formally, Ro — ► oo as 
cf) -> 0. 

In QG theory, adiabatic dynamics is governed by the material 
advection of potential vorticity: 



where D/Dt is the material derivative and 

f^^i ■ f2 d ( 1 d^\ 

is the QG potential vorticity in the '/3-plane approximation' (see 
section 2.3). Here / = fo + /3y, where fo = f((f>o) and f3 = 
{df /dy) 1 0=0 O for a specific latitude 0o; ip is the streamf unction; 
and, V 2 is the horizontal Laplacian operator. Note that q QG can be 
inverted - as with the full primitive equation q, under the QG bal- 
ance condition - to obtain all other dynamical variables (Hoskins, 
Mclntyre & Robertson 1985). 

The QG equations (the above advection equation for q QG plus 
boundary conditions) derive from the more complete primitive 
equations (e.g. Pedlosky 1987). The standard QG equations are 
appropriate for large-scale motions on many planets, away from 
the low latitudes. However, they are not broadly 1 appropriate for a 
large number of extrasolar planets, which are characterized by Ro 
of order unity (see Table 1) - even away from the equatorial region. 
More importantly, much dynamics of interest on extrasolar planets 
occur in the equatorial region (section 3.5), where the traditional 
QG theory does formally break down. Therefore, we perform our 
stability analysis using the full primitive equations. 



QG theory may still be valid locally on hot extrasolar planets. 
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2.2 Governing Equations 

In the standard pressure (p) coordinate system (Kasahara 1974), the 
hydrostatic primitive equations read: 



<9v <9v 

— + v-Vv + w— + /k x v + V$ J^v 
at 

80 ~ d6 

— +v-V0 + u; — = To 

at op 



V-v + 



duo 
dp 

+ h0 = 







(la) 
(lb) 
(lc) 
(Id) 



Here v(x, t) = (u,v) is the (zonal, meridional) velocity in the 
frame rotating with ft, where x £ M 3 ; a; = Dp/Dt is the 'ver- 
tical' velocity, where D/Dt = d/dt + v- V + cod /dp with V 
operating along constant surfaces of p (which in general is not ma- 
terially conserved); = T(p r /p) K is the potential temperature, 
where T is the temperature, p T is the reference pressure, k — 7Z/c p 
with 1Z the gas constant and c p the specific heat at constant pres- 
sure; $ = gz is the geopotential above the planetary radius R p , 
where g is the constant surface gravity; h{p) — 7l(p/p v ) K / p\ 
and, T w and Tq are, respectively, momentum and potential tem- 
perature sources/dissipations. From here on, we exclusively work 
in p-coordinate and drop the tilde over the advective and gradient 
operators for notational clarity. 

The above set of equations is closed with the ideal gas law, 
pa — VST. The equations are also supplemented with the boundary 
conditions, 



cj = 



at 



V = 5 Pr • 



(2) 



Hence, the domain boundaries are material surfaces and no mass 
flows across them. 



2.3 Two-Layer, Beta-Plane Analysis 

In this section, we linearize equations (1) on the /3-plane, where 
f(<p) is represented by f(y) = fo + fly. Here, fo and f3 are 



constants, y 



Rr> 



and the motion is assumed to be 



periodic in the zonal direction with no meridional component at 
the latitudinal boundaries. The /3-plane is a tangent plane located 
at 0o, and the setup is only formally justified for scales that are 
small compared to R p . However, in practice the /3 -plane approx- 
imation mainly results in small distortion of planetary waves and 
captures the essential qualitative behavior. For the analysis in this 
section, we neglect source/dissipation terms in equations (1) - i.e. 
T w (x, t) = Tq (x, t) = 0. This is because, as discussed in sec- 
tion 1, we are interested in the dynamics of jets that result from 
conditions e.g. when the net heating is not large or the effective 
thermal relaxation time is not small. 

We perform a standard normal mode analysis of the baroclinic 
instability admitted by a two-layer representation of equations (1). 
Similar work has been carried out by Wiin-Nielsen (1963) and 
Fraedrich & Frisius (2001) for the Earth. As in these studies, we 
simplify equation set (1) to that appropriate for a discretised model 
with two equally- spaced, stacked layers in the p-coordinate. In this 
model v, and $ are defined at odd levels and cj is defined at even 
levels. The structure is illustrated in Fig. 1. The equations for the 



Pressure Levels 
hPa — 



250 hPa 1 
500 hPa 2 



Variables 
■ co =0 



750 hPa 3 v 3 , e 3 , <£ 3 

1000 hPa 4 co 4 =0 

Figure 1. Vertical structure of the two-layer primitive equations model. 
Different field variables are defined on different levels; 1 hPa = 10 2 Pa 
= 1 mbar. Bold lines are the layer boundaries. 



interior levels are: 

dvi / v 3 - vi 

-¥ + VrVvi+W2 iw 

<9v 3 / v 3 - vi 

+ V 3 -VV 3 +CJ 2 

dt \ 2 Ap 

^— + V-(0i vi) + -— - 
at Ap 

d0 2, T-i ( r\ \ OJ2O2 ~ 



+ /kxvi = 
+ /kxv 3 = 



-V$i 

-v$ 3 



+ v-(0 3 v 3 ; 

V-vi + ^ = 
V-v 3 - ^ = 



dt v oy Ap 

Ap 
Ap 

$1 - $3 - h 2 Ap0 2 . 



(3a) 
(3b) 
(3c) 
(3d) 
(3e) 
(3f) 
(3g) 



where Ap — p r /2 denotes the pressure difference between odd or 
even numbered levels and 62 = = (0\ + 63) /2. It follows from 
equations (3e) and (3f) that barotropic (vertically averaged) wind 
is non-divergent. In the present analysis, we take the Brunt- Vaisala 
frequency N to be uniform; GCM simulations by Thrastarson & 
Cho (2010) show static stability to be fairly constant over one or 
two scale heights for a wide range of conditions. 

Baroclinic instability in the two-level primitive equations sys- 
tem is obtained from perturbations of an unstable 'Eady-type' basic 
flow with uniform vertical shear (Eady 1949): 



Ul 
Vi 



-u 3 - 

V 3 = 

0J3 = 

2/0 

h 2 Ap 
h 2 Ap 



U 



U y + cr 
U y - a 

Uoy. 



2/o 



h 2 Ap 

Here Uo (= U/2) characterizes the strength of the thermal wind 
and its shear; and, cfq = (61 — # 3 )/2 is related to the reference 
static stability, S = —a d hiQ/dp, through 

n cr /i 2 



For simplicity we shall consider meridionally-independent pertur- 
bations applied to the above basic flow. We now linearize equa- 
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tions (3) about this basic state and arrive at the following: 



— + U — 
dt dx 

dt dx 



Ap 

Uo_ 
Ap 



u 2 i + /k x vi - 
uj' 2 i + /k x v 3 - - 



dx 
Ox 



W foUo , . 1 v ^ 



= 



V-v 3 







Ap 

$i -$3 = h 2 Ap6' . 



(4a) 

(4b) 

(4c) 

(4d) 

(4e) 
(4f) 



The temperature equation (4c) is obtained by summing equa- 
tions (3c) and (3d) and linearising. Further, if we denote the vertical 
average of a variable £ by 

£ + = i(fc+&) 

and the half vertical difference by 

e-^«i-&), 

summing and differencing equations (4a) and (4b) give: 



dt dx 



Uo_ 
Ap 



uj' 2 i + /kx v+ 



<9x 



(5a) 



"aT + ^ "a^ + / kx v - = B— • ( 5b ) 

at ax ax 

By applying curl and divergence, we obtain vorticity and diver- 
gence forms, respectively, of the above equations. The equations 
set is closed when potential temperature and pressure velocity are 
eliminated using the hydrostatic and continuity equations. We then 
introduce the streamf unctions, ipi and ^3, and the velocity poten- 
tials, xi an d X3, for levels 1 and 3 such that 

d 2 

^(Xi+Xs) = 

and obtain four evolution equations for the barotropic vorticity, 
baroclinic vorticity, baroclinic divergence, and geopotential (i.e. 
potential temperature): 





= ^( 
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Ox 2 \ 
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respectively. The evolution equations for these quantities are: 
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dt 
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dx 



dx 2 



dx 



dx 2 dx 2 dx 

gV4 _ na d 2 X - 

dx dx 2 



(6a) 
(6b) 
(6c) 
(6d) 



At this point, the foregoing system of equations can be made 
non-dimensional for a more 'generalized' treatment, as is typical in 
instability studies. However, we shall describe our analysis of the 
equations presented in the dimensional form. We feel this facilitates 
a more lucid interpretation of the results in some ways. For the 
interested reader, we have included the non-dimensional account in 
Appendix A and refer the reader to that section, especially for the 
dependence of the results on non-dimensional parameters. 

Denoting disturbances by 



\l> = \frexp{i/c (x - ct)} , 



where* = (^,f ,x'_,$'_) T ,* 
c £ C, equations (6) reduce to 



and 



with 
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For a non-trivial solution, det(M) = 0. This leads to a fourth-order 
characteristic equation for c : 



+ c 



k 6 



30 
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Ml 

k 4 
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fSuS P 2 
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= 0. 



ul 



(7) 



Equation (7) is solved numerically for c as a function of k, 
while keeping the values of fo , f3 , Uo , 71 , (Jo and k constant. If 
3tti{c} / 0, the disturbances grow or decay exponentially since 
they are proportional to exp{— i kct}. Two of the roots of equa- 
tion (7) are stable eastward- and westward-traveling inertia-gravity 
waves. The other two roots are baroclinic waves. These waves prop- 
agate neutrally (i.e. without growing or decaying) eastward and 
westward, if Sm{c} = (provided 5Re{c} / 0). 

The baroclinic wave solutions to equation (7) are presented 
in Fig. 2 for a planet with HD209458b parameters given in Ta- 
ble 1. The left panel shows the growth rate, /c-$tti{c}, as a func- 
tion of wavelength 2nk~ 1 at several different latitudes: cj) — 
(60°, 45°, 35°, 25°). The growth rates are labelled 'HD60' (red), 
'HD45' (green), 'HD35' (yellow) and 'HD25' (blue), respectively. 
The right panel shows the corresponding phase wave speeds 5Re{c}. 
Note here that unstable baroclinic waves travel westward relative to 
the mean flow (9fte{c} < 0). 

From Fig. 2 it can be seen that the wavelength of the most 
unstable mode at = 60° is 1.7 x 10 8 m, corresponding to 
1.8 undulations (i.e. ~2 crests and troughs each) at this latitude. 
The growth rate of the instability is 3.1 r _1 , where r = 2tt Q _1 is 
the planetary rotation time. At cj) — 45° and cj) — 35°, the most un- 
stable modes correspond to 2.2 and 2.3 undulations around their re- 
spective latitude circles and with growth rates 2.3 r _1 and 1 .5 r _1 , 
respectively. Hence, jets centred at lower latitudes have increased 
growth times with modestly increased wavelengths of the most un- 
stable mode. Significantly, linear analysis predicts stability for jets 
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Growth Rate x 10 -6 s " 1 Phase Speed m s~ 1 




Wave-length x 10 8 m 



Figure 2. Growth rate [k • Sm{c}] (left) and phase speed [3£e{c}] (right) for HD209458b, as a function of wavelength 2ir fc -1 . Curves 'HD60', 'HD45', 
'HD35' and 'HD25' represent growth rates and phase speeds at <j> = (60°, 45°, 35°, 25°); f = 4.2 x 10 -5 sin</>s _1 ,/3 = 4.2 x 1CT 13 cos^m -1 s _1 , 
U = 500 m s _1 , n = 3500 J kg -1 K _1 , a = 300 K and k = 0.286. Curve 'HD45L' has been computed for HD209458b parameters at <j> = 45°, but 
with?7 = 200 ms -1 . 



Growth Rate x 10" 6 s " 1 Phase Speed m s -1 




0.2 0.4 0.6 0.8 1 



Wave-length x 10 8 m 



Figure 3. Growth rate [k • ^sm{c}] (left) and phase speed [3fte{c}] (right) at 4> = 45° for Earth, Jupiter and GJ436b as a function of wavelength 2tt k 1 . 
For the Earth, f = 10~ 4 s _1 , P = 1.6 X lO" 11 m -1 s -1 , U = 20 m s _1 , U = 287 J kg -1 K -1 , a = 15 K, K = 0.286. For Jupiter, 
fo = 2.5 x 10- 4 s- 1 , (3 = 3.5 x lO -12 m -1 s" 1 , U = 25 m s -1 , K = 3779 J kg -1 K _1 , ctq = 24 K, k = 0.286. For GJ436b, f = 3.9 x lO" 5 s" 1 , 
P = 1.4 x 10- 12 m- 1 s~\ U = 250 m s _1 , U = 3500 J kg" 1 K _1 , a = 150 K, k = 0.286. Note, a values for Jupiter and GJ436b have been 
computed assuming constant temperatures of 120 K and 750 K, respectively. Note the change in scales, compared with Fig. 2. 



located at or equator-ward of cj) = 28° (see e.g. the flat, blue curve 
labelled 'HD25'). 

To illustrate the dependence of the growth rate and phase 
speed on the characteristic flow speed (or, equivalently, shear 
strength), we also present in Fig. 2 result obtained for the case with 
U = 200 m s _1 at = 45° (black curve labelled 'HD45L'). Com- 
paring the 'HD45' (green) and 'HD45L' (black) curves, we see im- 
mediately that the growth rate of the most unstable mode decreases 
significantly for the smaller Uo case. The instability takes ~ 4 times 
longer to develop in the weaker speed/shear case. We also note 
that the wavelength of the most unstable mode decreases slightly. 



Hence, as Uo decreases, the number of undulations increases for a 
jet located at a given latitude. 

The qualitative behaviour described above is not restricted to 
HD209458b. It applies to any planet that has a meridional tem- 
perature gradient. To illustrate the general applicability of our re- 
sults, we present in Fig. 3 the growth rates and phase speeds at 
<j) — 45° for several planets: Earth, Jupiter and GJ436b (red, green 
and black curves, respectively). For the Earth, the wavelength of 
the most unstable mode is 4100 km, corresponding to ~ 7 un- 
dulations at midlatitude, with growth rate of 1.6 r _1 (i.e. growth 
time of 15 hours). This is consistent with many studies of baro- 
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clinic instability on the Earth (e.g. Thorncroft et al. 1993; Polvani et 
al. 2004). The corresponding values of undulations for Jupiter and 
GJ436b are ~ 43 and ~ 1 with growth rates 0.48 r" 1 and 0.56 t~\ 
respectively. Accordingly, if baroclinic instability occurs on these 
planets, it appears Jupiter simulations must be of very high resolu- 
tion to capture the instability. On the other hand, the instability at 
the midlatitude of GJ436b would clearly be of planetary scale and 
thus may lead to a possible observable variability signal for this 
planet on a timescale of ~ 1.8 planetary rotations. Note that the 
phase speeds of the unstable baroclinic waves on Earth and Jupiter 
are very small (close to zero) compared to those on the extrasolar 
planets, HD209458b and GJ436b. 

We have also carried out linear growth rate analysis with the 
two-layer QG model for HD209458b and have compared the results 
with those from the primitive equations model, presented above. 
In the two models, the growth rates at high latitudes and midlati- 
tudes are equivalent to within 5 per cent. However, at low latitudes, 
the QG model overestimates the growth rates by approximately 
25 per cent. Moreover, the QG model predicts instability down to 
4> — 23°, whereas the primitive equations model predicts insta- 
bility only down to (j) — 28°. Below these latitudes, both models 
predict stability. Thus, ageostrophy appears to provide a stabilizing 
factor in this case. Given that inertia-gravity waves are not filtered 
in the primitive equations model (as they are in the QG model), the 
enhanced stability may be due to the gravity waves 'leaking away' 
some of the energy that drive the instability. 

2.4 Limitations 

The preceding analysis is highly idealised. Therefore, it has lim- 
itations. For example, in general, planetary jets possess a three- 
dimensional structure - with concurrent vertical and meridional 
shears, as well as zonal asymmetry. Also, the atmosphere is con- 
tinuously stratified. One effect of a two-layer discretisation with 
uniform zonal flow in each layer is the inability to capture the sym- 
metry breaking between eastward and westward jets. These limita- 
tions are discussed more in detail below. 

In flows with both vertical and horizontal shears, the growth 
of unstable baroclinic waves may be suppressed by the 'barotropic 
governor' effect (e.g. James 1987; Nakamura 1993a; Pedlosky 
1964). The effect is not fundamentally related to the sign of the 
jet, but a key ingredient is a counter-gradient eddy momentum flux 
u'v' generated under a horizontally sheared flow; here the overbar 
indicates a zonal average. The shear and the momentum flux rein- 
force each other to distort the meridional structure of the wave, sup- 
pressing the growth rate and shortening the wavelength of the most 
unstable mode. Thus, the full non-linear evolution of the instabil- 
ity exhibits lower growth rates and shorter wavelengths, compared 
with those indicated by the linear analysis presented in this section 
(see section 3). 

The atmosphere is also continuously stratified. A representa- 
tion more realistic than a two-layer model changes the instabil- 
ity properties described in this section. The main change is that 
the short-wave and long-wave cut-offs in the two-layer represen- 
tation (see Figs. 2 and 3) no longer exist in the continuum of un- 
stable modes (e.g. Charney 1947; Green 1960; Kuo 1979). The re- 
tained modes (Charney and Green modes, discussed below) are not 
expected to change qualitatively the asymptotic behaviour of the 
instability. However, they do provide additional modes for wave- 
wave interaction during the non-linear growth phase - hence, affect 
the details of the evolution; this may be significant for finite-time 
variability. 



Another limitation of the linear model presented is that it does 
not distinguish between the signs of the jet (or shear). This is be- 
cause symmetry is preserved under the interchange of the shear 
sign, given the laterally uniform flow; hence, distinction between 
the two signs is not expected. This is in contrast to the flow used 
in the non-linear calculation (section 3), in which the growth rate 
for an unstable westward (negative shear) jet is smaller than that 
for the unstable eastward (positive shear) jet at the same latitude. 
The two signed flows behave differently in this case because of the 
change in the sign of the jet curvature. Furthermore, a westward jet 
has only one unstable mode (Charney mode) as opposed to an east- 
ward jet, which has an infinite number of unstable modes (Green 
modes). 

A similar observation has been made by Wang (1989), who 
observed a difference between eastward- and westward- sheared 
baroclinic flows in the Charney model (continuously stratified QG 
model on the /3-plane). He has pointed out that the maximum 
growth rate for the absolute value of non-dimensional shear is sub- 
stantially smaller for a flow with westward shear than a flow with 
eastward shear. Moreover, while the eastward jet is baroclinically 
unstable for any value of vertical shear A, the westward jet is un- 
stable only if 

Jo 

Note that, at = 45°, the critical shear for HD209458b parameters 
used in this work is:A c = — 7.9xl0 _4 s _1 . The shear of the unsta- 
ble westward jet described in section 3.4 is A = —1.7 x 10 -3 s _1 , 
consistent with (8). 



3 NON-LINEAR EVOLUTION 
3.1 Numerical Model 

To study the full non-linear evolution, we use a well-tested par- 
allel pseudospectral model, BOB 2 (Scott et al. 2003). This model 
solves equations (1) in spherical geometry, subject to the boundary 
conditions (2). As in many models, the equations in the vorticity- 
divergence form are solved, where (relative) vorticity is C(x, £) = 
k • V x v and divergence is £(x, t) = V • v. The equations in this 
form are more amenable for the spectral transform method (Orszag 
1970; Eliasen et al. 1970; Canuto et al. 1988). Domain decomposed 
spectral transform algorithm is used in the horizontal direction and 
standard second order finite difference scheme is used in the verti- 
cal direction. The latter direction is in pressure coordinates in the 
numerical model. 

To follow the evolution over long duration, explicit dissipation 
is applied so that artificial accumulation of energy at small scales is 
prevented (see e.g. Cho & Polvani 1996). The dissipation is in the 
form of a linear superviscosity operator, — z/V 4 (-), applied to the 
prognostic variables, {C, S, 0}; here v is constant. A small Robert- 
Asselin time filter e (Robert 1966; Asselin 1972) is applied, at ev- 
ery time step and in each layer, to filter the computational mode 
arising from using a second-order time-marching scheme (see e.g. 
Thrastarson & Cho 2011). No other numerical dissipators, fixers, 
stabilisers or filters are used in performing the simulations. 



2 'Built On Beowolf 
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3.2 Model Setup 

To study the non-linear evolution of baroclinic instability, we ini- 
tialize our model with an idealized jet that satisfies the necessary 
condition for baroclinic instability, the Charney-Stern-Pedlosky 
condition described in section 2. The jet is initially set to be ei- 
ther eastward or westward, and centred at a latitude between 0° 
to 60° N. A large number of simulations have been performed for 
this study, carefully varying each parameter (jet location, strength, 
shear, profile, direction as well as domain size, etc.) in an indepen- 
dent series of simulations. A very small subset of these runs, which 
we use for discussions in sections 3.2 to 3.5, is given in Table 2. 
The set illustrates the basic points we wish to make. 

All the jets are initially non-linearly balanced so that a self- 
consistent background temperature structure is generated (Fig. 4). 
The jets are then bumped at the beginning of the simulation by 
an infinitesimal temperature disturbance which is independent of 
altitude, a barotropic 'heat bump', and allowed to evolve freely 
thereafter. The setup is chosen to be similar to that in Polvani et 
al. (2004) for validation and comparison purposes. For example, 
following that work, the initial zonal flow uo in our runs is, in gen- 
eral, 



Here 



U sin m [7rsin 2 (0-0o)]F(^*), O < < fa 
, otherwise . 



1 - tanri 



• z 



Azq 



Zl 



(9) 
(10) 



with z* — —H\og(p/p r ), and 0o and 0t are taken to be the fol- 
lowing: 0o = and 0t = vr/2 for jets centred at midlatitude 
(E45N and E45N2b), O = vr/12 and T = tt/2 for jets cen- 
tred at 60°N (W60N) and O = -7r/4 and T = vr/4 for jets 
centred on the equator (EEQ). The typical values of the parame- 
ters are: U ±1000 m s~\ z 1823 km, zi 2486 km, 
Az = 414 km, H = 580 km, and p r = 10 5 Pa (= 1 bar). 
The latitudinal width of the jet is determined by m in (9), where 
m = 3 corresponds to a jet width of ~ 40° (Fig. 4a and Fig. 4b) 
and m = 1/2 to a width of ~ 85° (Fig. 4c). To discuss jets that 
closely match those produced in current GCM simulations of ex- 
trasolar giant planets, we present runs which are initialized with 
wider (m = 1/2) jets in the equatorial region and narrower jets 
(m = 3) poleward of 45° N. 

The basic state temperature, T = Tb(</>,p), is obtained by 
combining meridional momentum and hydrostatic balance equa- 
tions: 



dTo 



H 

a 



(R P f + 2uq tan c 



duo 



(11) 



Integrating (11) results in a temperature distribution that is in non- 
linear, gradient wind balance with the specified jet. Here we have 
used a reference temperature of 1500 K as the constant of integra- 
tion. The value is consistent with initial conditions and results of 
many GCM calculations. The basic state flow uo((j),p) and poten- 
tial temperature 0o (0, p) for runs E45N (eastward midlatitude jet), 
W60N (westward high latitude jet) and EEQ (wide eastward equa- 
torial jet) are shown in Fig. 4. Recall that 0q is related to T by 
0o = To(p r /p) K . To catalyse the instability, To is given a small 
perturbation T' in the form of a localized bump at all pressure lev- 
els such that 



Table 2. Summary of jet configurations discussed: m is a parameter that 
controls the jet width [see equation (9)]. Note, in run E45N2b the bot- 
tom boundary is set at p = 2 bar and the vertical structure function 
F(z*) in equation (9) is specified as = { 1 - tanh 8 [(z** - 

z 2 h) / Azq]} s\n A (n z** / z\) , where z** = -H log[(p+po)/Pr], with 
po = 60 hPa and z 2 b = 900 km. 



Run 



Width (m) Latitude Direction 



E45N 


3.0 


45° N 


East 


E45N2b 


3.0 


45° N 


East 


W60N 


3.0 


60° N 


West 


EEQ 


0.5 


0° 


East 



T'(A, 0) = .Asech 2 [3 (A - A )] seen 2 [6(0 - O )] , (12) perturbation. 



for — 7T < A < 7T. Here A = 1 K and (0o, Ao) represents the jet 
centre (latitude, longitude). 

Our vertical domain, which typically extends from 1 to 
10~ 3 bar, is resolved by 20 equally spaced pressure levels. The hor- 
izontal resolution of results presented in section 3.3 to section 3.5 is 
T170, or 170 sectoral modes and 170 total modes in the spectral ex- 
pansion (see e.g. Thrastarson & Cho 201 1). The resolution is des- 
ignated 'T170L20'. The inverse transformation is performed on to 
a 512x256 Gaussian grid covering the entire globe. The grid size is 
chosen for de-aliasing (Canuto et al. 1988, and references therein). 
Equations (1) are integrated for up to 60 r (i.e. 60 planetary rota- 
tions) with v — 6 x 10 19 m 4 s _1 . A timestep size, At = 30 s, 
and a Robert- As selin coefficient, e = 0.01, are used for the time 
integration. 

As already mentioned, the choice of our initial conditions is 
partly motivated by current GCM results of hot extrasolar giant 
planet atmospheres. These studies suggest typical flow speeds of 
O(100 — 3000 m s _1 ) and zonal flow consisting of up to ~ 3 jets - 
often a broad equatorial eastward jet and a smaller amplitude nar- 
rower westward jet at a higher latitude on both northern and south- 
ern hemispheres (e.g. Showman et al. 2008; Rauscher & Menou 
2010; Thrastarson & Cho 2010; Heng, Menou & Phillips 2011). 
The altitudinal and latitudinal profiles used here roughly mimic 
those presented in fig. 9 of Showman et al. (2008) and fig. 3 of 
Rauscher & Menou (2010). 

In what follows, we first describe the evolution of the midlati- 
tude eastward jet (run E45N). Although such a jet is not commonly 
observed in current simulations of hot extrasolar giant planets, re- 
viewing this case is useful because it allows the present work to be 
compared with analogous studies - and observations - of the Earth 
and because it allows a baseline to be constructed for other initial 
conditions presented here, namely the high-latitude westward and 
equatorial eastward jets that match more closely with aforemen- 
tioned extrasolar planet simulations. 



3.3 Paradigm Case 

Run E45N is the 'paradigm case'. It illustrates a typical non-linear 
evolution of a perturbed, marginally stable 3 baroclinic jet on a 



3 While the jets satisfy the (necessary) condition for instability, they re- 
quire an initial perturbation to evolve: they are perfectly stable without the 
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Figure 4. The basic state zonal wind uq [m s 1 ] (red) and potential temperature Oq [K] (black) as a function of latitude and pressure for runs (see Table 2): 
a) E45N, b) W60N, and c) EEQ. Contour interval for the zonal wind is 100 m s _1 and for the potential temperature 100 K. Negative contours are dashed. 



hot extrasolar giant planet in numerical simulations with high res- 
olution. The jet is zonally symmetric and eastward with speed 
1000 m s _1 at the jet core and decaying to zero at the periphery (see 
Fig. 4a). It meets the necessary conditions (i), (ii) and (iv) for baro- 
clinic instability, defined in section 2, with (i) and (ii) only weakly 
satisfied. This can be seen from Fig. 5, which shows (dqo/dy)e 
evaluated on an isentrope as a function of 4> and p. Note, here go is 
the potential vorticity defined on isobars, 

go(<M = -p(/k + Vxv )-V(9o , 

where V is the three-dimensional gradient operator in (A, <j>,p) 
space; and, (dqo / dy)e is a derivative taken along an isentrope such 
that 

\dy) e \dy) p \ dy ) p \ dp J dp' 

where y = R v (j) and [d( • )/dy] p is the derivative taken on an iso- 
bar. Other cases, with jets of different sign or location, are to be 
compared with this one. Fig. 6 presents the evolution of T (left col- 
umn) and C (right column) fields at the p = 975 hPa surface from 
run E45N, for r — to r — 8. The fields near the reference pres- 
sure level are shown since the kinetic energy is the maximum at 
the lower boundary for jet profiles shown in Fig. 4, similar to Gall 
(1976) and Simons (1972). Note that for these jets T « at this 
pressure level. 

In Fig. 6, the perturbed jet undergoes initially a period of lin- 
ear growth (r <C 4), when the most unstable mode emerges. At 
this early stage, the T field shows a small- amplitude perturbation 
from zonal symmetry. The £ field, on the other hand, is much more 
dynamic. At r — 4, finite- amplitude wave breaking in the £ field is 
already clearly evident, and the perturbation in this field is charac- 
terized by a distinct northwest-southeast tilt on the poleward side 
of the jet and southwest-northeast tilt on the equatorward side of 
the jet. The enhancement of the tilt proceeds concomitantly with the 
barotropic component of the flow, which generates negative merid- 
ional flux of the eddy zonal momentum (i.e. u'v' < 0) on the pole- 
ward flank of the jet and positive meridional flux of the eddy zonal 
momentum (i.e. u'v' > 0) on the equatorward flank of the jet (see 
e.g. Nakamura 1993b). 

During this early stage of the evolution, conversion of avail- 
able potential energy (APE) into eddy kinetic energy (EKE) slowly 




30N 



Figure 5. Meridional cross-section of the meridional potential vorticity 
gradient (dqo/dy)o for run E45N (northern hemisphere). Maximum and 
minimum values are ±3 x 10 -12 Km kg -1 s _1 with contour interval 
2 x 10 — 13 Kmkg - 1 s — 1 . Negative values are in blue and positive are in 
red. The zero contour is drawn with double thickness. 



begins, as can be seen from Fig. 7. These quantities are defined: 

EKE = [Wf + (v) 2 } dp dA , (14) 

where 0" is the deviation of from its isobaric average 0, and the 
integrations are over the surface area A and pressure p. The baro- 
clinic instability taps the APE to drive the eddy motions. 

At r « 5 a rapid non-linear development ensues in both fields. 
Note, for example, the scale change in the ( field at r = 6. A large 
amplitude wave can now also be clearly seen in the T field. In both 
fields, sharp frontal features form. These dynamically-significant 
sharp features require very high resolution to capture faithfully. 
This will be discussed more in detail in section 3.6. By r = 8 
sharp temperature gradients trail out around the anticyclonic re- 
gion (large 'clover-leaf shaped area of negative vorticity, shaded 
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Figure 6. Temperature T (left) and relative vorticity £ (right) from run E45N in polar stereographic view, centred on the north pole. The fields are shown at 
the 975 hPa pressure level for r = to r = 8. Maximum and minimum values for T are 1280 K and 1520 K, respectively, with contour interval 6 K. For f , 
the maximum and minimum values are ±5 x 10 -7 s _1 (r = 0), ±1 x 10 -6 s _1 (r = 4), ±1 x 10 -5 s _1 (r = 6) and ±4 x 10 -5 s _1 (r = 8); the 
contour intervals are, respectively, 2x 10 -8 s -1 ,4x 10 _8 s _1 ,4x 10 _7 s _1 and 1.6 x 10 -6 s -1 . The spectral resolution of this simulation is T170L20 
(see text). Note the large, nearly two orders of magnitude, change in the magnitude of £ during the evolution - as well as the formation of sharp fronts and 
coherent vortices, particularly at r = 6 and r = 8. 
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Figure 7. Evolution of globally-averaged eddy kinetic energy (per area) 
[J m~ 2 ] for runs E45N (dashed line), W60N (dotted line) and EEQ (solid 
line). The eddy kinetic energy for EEQ has been multiplied by a factor of 
50. 



in blue), forming curved baroclinic fronts. At this time, the EKE 
is well into its non-linear growth stage. Note the pools of warm air 
that have been pinched off (cyclonic vortices embedded in the anti- 
cyclonic region), intruding into the high latitudes. Simultaneously, 
broad regions of cool air spread into the tropical region from higher 
latitude (i.e. negative heat transport). Thus, the original equator- 
pole temperature gradient is significantly reduced by the instability. 

The poleward heat transport can be checked against linear the- 
ory for eddy transport (see e.g. Holton 1992; Vallis 2006) by ex- 
amining the zonal cross-sections of streamfunction, meridional ve- 
locity and temperature perturbations: ip\ v and T' ', respectively. 
The cross-sections at midlatitude are shown in Fig. 8. As can be 
seen, ip' and v tilt westward with height and T' tilts eastward with 
height, demonstrating that heat transport is taking place. Note that, 
in the case of baroclinically unstable westward jet at the same lati- 
tude, the directions of the tilts are reversed. This is because gradient 
wind balance produces a temperature distribution that is warmer at 
the poles than at the equator, as can be seen in Fig. 4b. This results 
in an equatorward transport of heat by the eddies (section 3.4). 

The long-time evolution of the run presented in Fig. 6 is illus- 
trated in Fig. 9 (r = 10 and r = 40). By r « 40 the T and C 
fields have organised into essentially zonal structures and eddy ac- 
tivity has mostly ceased. The cyclones that have emerged from the 
baroclinic wave breaking, strongly interact (r = 10) and ultimately 
merge into an unsteady cyclonic polar vortex (r = 40). A similar 
'end-state', resulting from vortex mergers, has been observed in 
HD209458b simulations of Cho et al. (2003). 

The long-range interaction of the like- signed vortices on hot 
extrasolar planets is more pronounced than on the Earth (and other 
cool, rapidly-rotating planets). This can be explained by the much 
larger Rossby deformation length scale, L^/R p = (9(1), on the 
hot extrasolar giant planet. Larger Ld means more robust mergers 
and a more dynamic final vortex (Cho et al. 2003, 2008; Cho & 
Polvani 1996). Scott (2011) has recently quantified this behaviour: 
merger and poleward migration of cyclones ensues if the potential 
vorticity anomaly q associated with a vortex exceeds the magni- 
tude of the planetary vorticity 2Q by ~ 12 per cent. In our case, 
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Figure 8. Perturbation streamfunction ip' (top), perturbation meridional ve- 
locity v' (middle) and perturbation temperature T' (bottom) at midlatitude 
as a function of pressure and longitude at r = 6 for run E45N. Contour in- 
tervals are: -24 x 10 8 m 2 s -1 to 11 x 10 8 m 2 s -1 in steps of 10 8 m 2 s -1 , 
-180 m s -1 to 320 m s -1 in steps of 20 m s -1 and -75 K to 125 K in 
steps of 5 K, respectively. Note, ip f and v' tilt westward with height and 
T' tilts eastward with height, signifying meridional transport of heat and 
reduction of equator-pole temperature gradient. 



we find (q' /2Q) ^ 1.19 - i.e. anomaly excess of 19 per cent - by 
r = 6, consistent with Scott's finding. 

The temporal evolution of the global average EKE (the dashed 
line for run E45N and solid line for run EEQ in Fig. 7) is typi- 
cally described as a 'baroclinic growth - barotropic decay' cycle 
(e.g. Simmons & Hoskins 1979; Thorncroft et al. 1993). During 
the cycle, conversion of APE to EKE is impeded by a positive 
feedback between the horizontal shear in the flow and the eddy 
momentum flux. At r ^ 10, the disturbances in run E45N are 
sheared out and EKE is lost to the mean flow through the Reynolds 
stresses more than it is gained through baroclinic conversion. The 
feedback is the main component in the previously mentioned non- 
linear 'barotropic governor effect', affected by the horizontal shear 
in the jet, spherical geometry and ageostrophy (see e.g. Nakamura 
1993b). 

The zonal mean zonal wind u and zonal mean zonal potential 
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Figure 9. Same as Fig. 6 but for r = 10 and r = 40. Contour interval for temperature is 6 K. The maximum and minimum contours for relative vorticity are 
±4 x 10 -5 s _1 at r = 10 and ±10 -5 s _1 at r = 40. The respective contour intervals are 1.6 x 10 -6 s _1 and 4 x 10 -7 s _1 . 
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Figure 10. Zonal mean zonal wind u (red) and potential temperature (black) contours for runs a) E45N at r = 30, b) W60N at r = 60 and c) EEQ at 
r = 60. Wind contour interval for runs E45N and EEQ is 100 m s _1 and for run W60N is 50 m s 1 . Temperature contour interval is 100 K. Negative 
(westward) wind contours are dashed. 
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temperature at the end of the life-cycle is presented in Fig. 10 
(see panel a on the left for run E45N). The jet itself has becomes 
broader and more barotropic - much like the 'LCI life cycle' re- 
ported in Thorncroft et al. (1993). Furthermore, the meridional en- 
tropy gradient d#/d(/> is significantly reduced compared to the ini- 
tial state, particularly at the lower levels in the domain (cf. Fig. 4a). 
Much of the APE is taken up by the kinetic energy of the zonal 
mean flow and the flow is accelerated there. To quantify the accel- 
erations, consider the transformed Eulerian-mean zonal momentum 
equation (e.g. Andrews & Mclntyre 1978): 



du 
~dt 



1 



d 



R v cos 6 d c 



(u cos <fi) — f 



du _* 
dp 

- V F. 



Here v* and oj* represent the 'residual' mean meridional circula- 
tion, 



_ _ d_ ( v'6' \ 
V ~ dp [dS/dp) 



1 d ( v'6' \ 
R p cos cj)d(j) \dO/dp COS(t) ) 



and F = (F</>, F p ) is the Eliassen-Palm (E-P) flux vector with 



— / v'O' \f du 

-U'V + -7T- 



Fp — Rp ( 



(C + /) 



\dO/dp)\dp 
dO/dp J 



The influence of eddies on the mean flow is measured by the 
E-P fluxes: a convergent flux (V • F < 0) corresponds to the de- 
celeration of the eastward flow and a divergent flux ( V • F > 0) 
corresponds to acceleration. Fig. 11 depicts vertically and tempo- 
rally averaged E-P flux divergence for E45N (dashed line) over the 
life-cycle. The E-P fluxes are divergent in the net on the poleward 
flank of the jet, where the flow is accelerated, and (more strongly) 
convergent in the net on the equatorial flank, where overall the flow 
speed is reduced from the initial value (cf. Fig. 10a ). 

Finally, we note that the most unstable mode for this calcu- 
lation is ~ 4 (see Fig. 6). As discussed earlier, the linear theory 
of section 2.3 underestimates this number to ~ 2. However, the full 
numerical simulation shows that the simple linear theory is success- 
ful, at least qualitatively, in capturing the behavior of the instability 
in the following sense: the most unstable mode and the growth time 
for the baroclinic wave amplitude for HD209458b are smaller than 
the corresponding quantities for the Earth (cf., for example, Polvani 
et al. 2004). 



3.3. 1 Lower Boundary 

As is well-known, boundary conditions are crucial in solving dif- 
ferential equations. Differences in the conditions, even in relatively 
simple physical situations, can alter the admitted solutions. For ex- 
ample, new or modified modes can be introduced or existing modes 
can be filtered by employing rigid boundary condition (i.e. w = 0). 
The lower boundary of the simulations discussed in this paper is 
rigid and located at 1 bar for the most part. In this case, the ver- 
tical wind shear in the basic flow used is small, but non-zero, at 
the bottom boundary and baroclinically unstable modes can arise 
due to the presence of the boundary - via condition (ii) of the 



0) 

o 

0) 
tuD 

?H 

CD 
! 



~i 1 1 r 

--W60N 

— E45N 
— EEQ 




-3 

X 
I 

Oh -6 

H -90 -60 -30 30 60 90 

latitude 

Figure 11. Vertically and temporally averaged divergence of Eliassen-Palm 
(E-P) flux [m s~ 2 ] for runs E45N (dashed line), W60N (dotted line) and 
EEQ (solid line) during the life-cycle in each run. The EEQ curve has been 
multiplied by a factor of 100. 



Charney-Stern-Pedlosky criteria. However, the stably- stratified ra- 
diative zone in hot extrasolar planet atmospheres may extend down 
to perhaps as deep as 1000 bars (Guillot & Showman 2002). Hence, 
the effects of lowering the bottom boundary and preventing flow 
shear there require careful consideration. 

Fig. 12 presents a run (E45N2b) that is very similar to the 
'paradigm case', but with the lower boundary of the calculation 
extended down to 2 bars. The jet is confined to pressures above 
the 1 bar level. In doing so we remove the influence of the lower 
boundary far enough away from the jet while still retaining an ad- 
equate vertical resolution. In the figure, uo and Oo are shown in 
the left panel. Note, the jet profile shown in the figure has a differ- 
ent vertical structure than the 'paradigm case' jet. This is because 
balancing a jet with vertical structure given by equation (9) to the 
isothermal reference state produces a convectively unstable region 
in the computational domain, causing the vertical coordinate to lose 
single- valuedness and the run to immediately crash. The shown 
profile does not suffer from this. Significantly, the Charney-Stern- 
Pedlosky criteria (iv) and (i) remain satisfied for this profile. It is 
crucial to understand here that, once these criteria are met, it does 
not matter whether the lower boundary is located at 10 or 1000 bars 
for the instability to occur. 

The C field at r = 10 is shown in the right panel of Fig. 12. 
Although the evolution is now slightly altered from the 'paradigm 
case' (i.e. mode-3 is dominant, rather than mode-4), the jet is still 
unstable, as expected. We have verified that the evolution in run 
E45N2b is indeed a result of baroclinic instability: the perturbation 
fields tilt in the appropriate directions with height, as seen in Fig. 8 
for run E45N. The instability is, however, weaker and evolves dif- 
ferently than when there is an initial vertical wind shear and merid- 
ional entropy gradient at the lower boundary. We have performed 
a simulation with the initial flow profile used in run E45N2b but 
with bottom raised to 1 bar, in which the shear and gradient is non- 
zero at the bottom. The peak global eddy kinetic energy in this run 
is ~ 40 per cent greater and vorticity perturbations are up to six 
times stronger than in run E45N2b. Nevertheless, the point is, the 
instability is present regardless of vertical flow shear at the bottom 
boundary. 
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Figure 12. Left: Basic state zonal wind uq (red) and potential temperature 0q (black) for run E45N2b with the lower boundary extended to 2 bars. Contour 
interval for the zonal wind is 100 m s _1 and for the potential temperature 100 K. Right: Relative vorticity £ from run E45N2b at r = 10. The field is plotted 
at 975 hPa pressure level. Maximum and minimum values are ± 6 x 10 -6 s _1 , with contour intervals of 4 x 10 -7 s _1 . 



3.4 High Latitude Westward Jet 

Having presented the evolution of the 'paradigm case', we now 
present the case of baroclinically unstable high-speed westward jet 
at high latitude (run W60N). The speed at the core of the jet is 
1000 m s _1 . Such jets have been observed in recent GCM simula- 
tions (e.g. Showman et al. 2008; Menou & Rauscher 2009; Thras- 
tarson & Cho 2010; Heng, Menou & Phillips 2011). In these sim- 
ulations, the high latitude jets tend to be more narrow and shallow 
than the equatorial jets. Equatorial jets will be discussed in sec- 
tion 3.5. The T field from run W60N at r = 9 (975 hPa pres- 
sure level) is shown in the left column of Fig. 13. Polar stereo- 
graphic view, centred on the north pole (top frame), and cylindrical- 
equidistant view, centred on the equator (bottom frame), are shown 
for latitudes poleward of 20° N. For comparison, the right column 
shows the corresponding projections of the T field from the E45N 
run at t — 7, roughly at a similar stage of the evolution in run 
W60N. 

In run W60N, baroclinic wave develops a significant 
northwest-southeast tilt. This is consistent with the predominantly 
negative momentum fluxes on the poleward side of the jet during 
the linear stage of the evolution. Again, the flow is characterized by 
sharp cyclonic fronts, this time with the most unstable mode having 
3 undulations at (fi = 60°. The reduction in the number of undula- 
tions is also consistent with linear theory developed in section 2.3. 
However, contrary to the predictions from the linear analysis, the 
growth rate of the instability is lower than for the E45N run (cf. the 
onset of growth between dashed and dotted lines in Fig. 7). As al- 
ready discussed, this agrees with the analysis of the Charney model 
by Wang (1989) and extends that result to the more general, global 
primitive equations model. 

Note that wave-breaking in the westward jet case occurs in 
the opposite direction to that in the 'paradigm case' (see bottom 
row of Fig. 13). The waves in run W60N (left) breaks eastward, 
whereas in run E45N (right) the waves break westward. 'Blobs' 
of higher temperature fluid penetrate into the lower temperature re- 



gion and cooler fluid subsides into the warm region. The situation is 
analogous to Rayleigh-Taylor or convective instability (e.g. Sharp 
1984), where a decrease of potential energy results under the inter- 
change of two blobs at different heights. In baroclinic instability, 
this can occur despite the stable density stratification because the 
density surfaces are sloping more steeply than the line joining the 
two blobs. Indeed, for this reason baroclinic instability is some- 
times refer to as 'sloping convection' (e.g. Vallis 2006). 

Note that the EKE for the westward jet case does not follow 
a simple 'baroclinic growth - barotropic decay' cycle, as seen in 
the 'paradigm case' (Fig. 7). Instead, after the initial decay stage 
at t ~ 15, the EKE for run W60N shows large vacillations, corre- 
sponding to a sequence of baroclinic-barotropic life cycles. Similar 
behavior of energetics has been observed by Feldstein (1991) for 
the Earth case. In a two-layer QG /3-plane model, he found west- 
ward jets to undergo a series of mixed, baroclinic-barotropic insta- 
bility, caused by the reversal of sign in the jet curvature d 2 m/dy 2 . 
Recall that /3 ^ 0. Hence, a barotropically unstable region, in 
which p — d 2 m / dy 2 < 0, forms at the core of the westward jet (as 
is the case in run W60N). The combined effects of vertical and hor- 
izontal shears reinforce each other to establish a mixed, baroclinic- 
barotropic unstable region. According to WKB analysis (e.g. Ben- 
der & Orszag 1999), growing disturbances emanating from a west- 
ward jet are trapped (i.e. reflected) between two turning latitudes, 
initiating the sequence. Consistent with this, the meridional struc- 
ture of the disturbance is able to remain close to the normal mode 
form. In contrast, disturbances emanating from the eastward jet are 
absorbed at or near the critical latitudes, resulting in a single cycle 
and meridional structure that changes with time. 

Fig. 10b shows the equilibrated u and at the end of the sim- 
ulation. The original westward jet has been completely disrupted, 
giving way to a fairly barotropic eastward jet centred at 60°. Pre- 
dominantly westward flow is now situated in the subtropics, at the 
upper levels of atmosphere. The reversal of the flow direction is 
consistent with E-P flux divergence shown in Fig. 1 1 (dotted line), 
which acts as a positive momentum source. Given that high-latitude 
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Figure 13. Temperature field T for run W60N at r = 9 (left) and T for run E45N at r = 7 (right). The fields are shown at 975 hPa pressure level. Top frame 
shows the field in polar stereographic view, centred on the north pole, and bottom frame shows the field in cylindrical-equidistant view, centred on the equator. 
In all frames area poleward of <f> — 20° is shown. For run W60N the maximum and minimum values are 1500 K and 1780 K, respectively. For run E45N the 
maximum and minimum values are 1280 K and 1520 K, respectively. Contour interval is 6 K in both runs. 



westward jets appear to be a fairly common feature in GCM simu- 
lations of hot giant extrasolar planets, the result here suggests ex- 
ternal (e.g. stellar irradiation) or internal (e.g. wave) forcing may 
be required to maintain baroclinic westward jets. Note also from 
Figs. 10b and 1 1, the negative zonal flow and E-P flux convergence, 
especially in the equatorial region. Significantly, such negative E-P 
flux divergences present a source of drag for equatorial jets. Finally, 
as in the 'paradigm case', the potential vorticity anomalies exceed 
2Q by over 12 per cent and by r — 35 a cyclonic polar vortex 
forms that is warmer than its surroundings (not shown). 



3.5 Equatorial Jet 

The evolution of a broad, high-speed equatorial jet (run EEQ) is 
presented in this section. The initial flow and potential tempera- 
ture is shown in Fig. 4c. Unlike the jets discussed in sections 3.3 
and 3.4, the equatorial jet satisfies the Charney-Stern-Pedlosky in- 
stability criteria (iv) on its flanks (at ~ 30° on both northern and 
southern hemispheres), rather than at the core. The stability of this 
jet's core is consistent with linear theory of section 2, which pre- 
dicts no growth for a jet located equatorward of 28°. 

Fig. 14 shows temperature T and relative vorticity ( fields at 
r — 26. At this time the instability is well developed, with sharp 
fronts rolling up non-linearly into cyclones at « 35°, where the 
instability criteria is met. A mode with ~ 7 undulations can clearly 



be seen at this stage of the evolution. The number of undulations 
is significantly higher and the growth rate is significantly lower for 
EEQ than for simulations where the same jet is placed at cj) = 30°. 
Evidently, since the vertical shear of the equatorial jet at its flanks 
is significantly lower than at its core, a value smaller than the peak 
core value for the characteristic flow speed should be used in inter- 
preting the results from the linear analysis. We have already seen 
that a weaker jet (shear) results in a smaller growth rate and wave- 
length of the most unstable mode at a given latitude (e.g. curves 
'HD45L' and 'HD45' in Fig. 2). Hence, our non-linear calculations 
appear to be in qualitative agreement with linear theory. Despite the 
instability at the flanks, the core of the jet in run EEQ remain stable 
throughout the integration (up to r = 60), in very good agreement 
with linear theory. 

The EKE evolution for run EEQ is shown in Fig. 7 (solid line). 
The equatorial jet instability is shallow and confined to a pressure 
range between 1 to 0.7 bar, unlike the high-latitude jet instability; 
in those cases, the range of instability is much larger, extending up 
to 0.01 bar. Thus, only the lower pressure levels exhibit an increase 
in EKE during the linear stage. For this reason, the EKE values 
have been multiplied by a factor of 50 in the figure: the globally 
averaged EKE for run EEQ is much lower than for the 'paradigm 
case' or run W60N. Qualitatively, the non-linear evolution of run 
EEQ is much like that of run E45N, with waves tilting and break- 
ing in same directions. However, potential vorticity anomaly only 
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Figure 14. Temperature T (left) and relative vorticity £ (right) for run EEQ in polar stereographic view, centred on the north pole. The fields are shown at 
975 hPa pressure level for r = 26. Maximum and minimum values for temperature are 1460 K and 1600 K, respectively, with contour interval 3.5 K. Values 
for £ are in the range ±1.5 x 10 _ 5 s -1 , with contour interval 6 x 10 -7 s -1 . 



slightly exceeds the polar value and the cyclonic drift does not en- 
sue. Therefore, a monolithic polar vortex does not form in this case. 

Interestingly, the jet structure is only slightly altered by baro- 
clinic instability from the basic state zonal flow. Mainly, the jet has 
become more barotropic at the flanks. This can be seen in Fig. 10c, 
which shows u and 6 at r = 60. Relatively small values of E-P flux 
convergence equatorward of 30° (see the solid line in Fig. 11) do 
not significantly contribute to the deceleration of the zonal mean 
zonal wind. 

It is also worth noting that we have investigated stability prop- 
erties of the westward equatorial jet and found it to be stable to 
baroclinic instability, in good agreement with Wang (1989). A 
westward jet placed at the equator would have to exceed sound 
speed, if the condition (8) of section 2 is to be fulfilled. However, 
we note that a broad, 'supersonic' westward jet does not appear to 
be unstable in full, non-linear GCM calculations (Thrastarson, pri- 
vate communication). But, the calculation is at T21 resolution (see 
next section). We have found eastward equatorial jets to be stable 
to baroclinic instability, if the width of the jet is 50° and smaller. 

3.6 Numerical Convergence 

Baroclinic instability in numerical simulations of hot extrasolar 
planets is highly sensitive to numerical resolution (both horizon- 
tal and vertical) and to dissipation. High resolution is required to 
capture the instability accurately. In particular, for the jet profiles 
used, five or more layers is necessary to capture the instability, and 
good convergence is reached only with ~ 10 or more layers. In ad- 
dition, high horizontal resolution ( > T85) is necessary to ensure 
accurate representation of the eddy fluxes, as well as convergence. 
Separately, artificial viscosity must not be too high, as it results in 
artificially-enhanced stabilization of the baroclinic modes. We em- 
phasise that resolution and dissipation requirements are dependent 
on the jet profile. Hence, the requirements should be carefully as- 
sessed for each profile employed. This 'problem-dependence' con- 



clusion has also been discussed by Thrastarson & Cho (2011) for 
'spin-up' experiments of hot giant planet circulations. 

Before presenting the results, a brief discussion concerning 
our general approach to convergence testing is in order. In gen- 
eral, for numerical stability reasons, the usual practice is to use a 
larger dissipation coefficient value when performing a calculation 
with lower resolution - or, alternatively, a smaller coefficient value 
when performing the same calculation at higher resolution - so that 
the damping time is same at the truncation scale. This results in a 
different damping time for a given mode at different resolutions. 
However, here our aim is to demonstrate convergence of the nu- 
merical model. Hence, we employ the same value of the coefficient 
at all resolutions so that each mode, up to the truncation, expe- 
riences same dissipation rate in all the runs. The employed value 
is: v — 6 x 10 19 m 4 s -1 . Similar methodology has been imple- 
mented in e.g. Polvani et al. (2004) to test convergence in the Earth 
case. Later, we also demonstrate non-convergence when the damp- 
ing time is chosen such that it is same at the truncation scale for all 
resolutions, as in the usual practice. 

The requirement of adequate resolution is demonstrated in 
Fig. 15. The figure shows a set of four simulations with all parame- 
ters identical to run EEQ presented in section 3.5 - except the reso- 
lution. The resolutions are: T21 (a), T42 (b), T85 (c) and T170 (d). 
Note, panel (d) is run EEQ. All four runs use the same value of 
superviscosity coefficient, v — 6 x 10 19 m 4 s _1 , as already men- 
tioned. Note also that the resolutions correspond, respectively, to 
64x32, 128x64, 256x 128 and 512x256 Gaussian grids in physi- 
cal space. But, because of the exponential convergence property of 
the spectral method, they are equivalent in accuracy to finite dif- 
ference grids O(10) times finer in resolution (e.g. Thrastarson & 
Cho 2011). Polar stereographic projections of the relative vorticity 
field £ at r — 22 are shown, when the instability is in the early 
exponentially-growing stage (see Fig. 7). 

Visual inspection of the fields readily reveals that the T21 (a) 
and T42 (b) runs do not converge to the T170 (d) run. The T85 (c) 
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Figure 15. Polar stereographic view of the relative vorticity field £, centred on the north pole, for four runs with all parameters identical - except the horizontal 
resolution. The common parameters are as in run EEQ. The fields at r = 22 are shown. The number at upper right in each panel indicates the resolution. 
Contour levels are the same to those in Fig. 14. 
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Figure 16. Kinetic energy density [m 2 s 2 (per wave number)] as a function of (total) wavenumber n. Spectra for the fields from the four runs shown in 



Fig. 15 (left). The different lines refer to different horizontal resolutions, as indicated in the legend. The viscosity coefficient is same (y = 6 x 10 9 m 4 s" 



in all four runs. Spectra for run EEQ at T 170 resolution with different viscosity coefficients (right): v - 
(black line), respectively. 



: 6 x 10 19 m 4 s _1 (red line) and v = 10 21 m 4 s 
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Figure 17. Same as Fig. 15 but with r d 1 = 5.07 x 10 4 s 1 in all runs. Maximum and minimum values for frames a) and b) are ±3 x 10 6 s 1 and for 
frames c) and d) ±1.5 x 10 -5 s _1 ; the contour intervals are, respectively, 1.4 x 10 -7 s _1 and 6 x 10 -7 s _1 . 



run is marginally converged, though this may change after a long 
time (e.g. many hundreds of planetary rotations). In the figure, 
frames (a) and (b) are qualitatively different than frames (c) and (d), 
which clearly show mode-6 instability. The T85 run in frame (c) 
captures the basic structure present in the T170 run in frame (d). 
However, spurious small-scale oscillations are also clearly visible 
in frame (c); these are not present in frame (d). The small-scale 
oscillations contaminate the calculation - causing the calculation 
to blow up, depending on the numerical parameters used; see e.g. 
discussion in Thrastarson & Cho (201 1). 

The above behavior can be quantified by computing the cor- 
responding kinetic energy spectra for each run. The spectra for the 
fields shown in Fig. 15 are presented in the left panel of Fig. 16. 
Inspection of the T85 and T170 spectra (red and black lines, re- 
spectively) confirms the convergence of the simulations. Note the 
presence of a clear dissipation range in the T170 run. In contrast, 
the appearance of nearly grid-scale waves in physical space for T21 
and T42 resolution runs corroborates the tendency of the spectrum 
(blue and pink lines for T21 and T42, respectively) in these runs to 
peel off and curl up considerably left of the aliasing limit (~ 21 for 
T21 and ~ 42 for T42). This is caused by discretization errors that 
are not adequately controlled by the applied explicit viscosity. 

We have also performed an analogous series of runs in which 
a much larger dissipation coefficient value, v — 10 21 m 4 s _1 
has been used. This mimics 'properly' dissipated runs at T21 and 
T42 resolutions (i.e. runs with well represented dissipation range 



in spectral space). However, in this series the high resolution calcu- 
lations are significantly over-dissipated and the physical space pic- 
ture is characterised by a severe reduction in eddy kinetic energy at 
all times. This is supported by the spectra for two T170 resolution 
simulations with the two coefficients (right frame of Fig. 16). The 
spectrum for a run using v — 10 21 m 4 s _1 (black line) is shown to- 
gether with the spectrum of the previously presented T170 run with 
^ = 6x 10 19 m 4 s _1 (red line). With the larger z/, the spectrum is 
severely over-dissipated with only ~ 20 modes being resolved; the 
rest of the modes clearly lie in the dissipation range. In contrast, at 
least rsj 80 modes are well-resolved with the smaller v value. 

It is important to understand that wavenumbers short- ward of 
the fiducial 'dissipation range' (i.e. less than 20 and 80, respec- 
tively, in the runs discussed above) are still affected by a small 
amount of dissipation in practice: that is, dissipation affects the 
entire spectrum of wavenumbers continuously, rather than just the 
wavenumbers in the dissipation range. The amount, while small, 
can nevertheless be dynamically significant, as it can change the 
quantitative character of the instability, even suppress the insta- 
bility altogether. Indeed, if the value of v is increased further, to 
10 22 m 4 s _1 , the baroclinic waves completely disappear. 

As discussed above, a common practice in numerical studies 
which vary the resolution is to adjust the dissipation coefficient v 
so that the damping time Td is same at the truncation scale n t : 

i 2 



Td 



n t (n t + 1) 



(15) 
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While the physical basis of this procedure is arguable, we demon- 
strate here that the practice still does not lead to convergence at low 
resolution. Consider the value, v — 6 x 10 19 m 4 s _1 , used in the 
high resolution simulation discussed above. The dissipation rate at 
the truncation scale for HD209458b corresponding to this v value 
isr^ 1 = 5.07xl0 -4 s -1 . Note that this damping rate is compara- 
ble to the rates used in current flow modeling studies of hot gaseous 
extrasolar planet atmospheres at resolutions lower than T170 (e.g. 
Rauscher & Menou 2010; Thrastarson & Cho 2010, 2011; Heng, 
Menou & Phillips 201 1). Its significance can be seen in Fig. 17. 

Fig. 17 shows the effects of adjusting v so that the damping 
rate is constant at the truncation scale. The rate used is the one 
just discussed above: r^ 1 = 5.07 x 10 -4 s _1 . Relative vorticity 
field C at four resolutions for run EEQ at r = 22 is shown. Two 
points are clear from the figure. First, the simulations in this series 
are not converged. The T21 and T42 resolution runs are completely 
over-dissipated and the momentum and heat transports are absent 
throughout the duration of the runs, up to r = 60. Second, the v 
values used in current GCM modelling studies of hot extrasolar gi- 
ant planets do not permit the instability. Note that, if the dissipation 
rate is chosen instead to be the one that 'adequately' permits the 
instability in the low resolution run, the high resolution runs are 
severely under-dissipated and inundated with noise (not shown). 
Either way, convergence is not achieved by fixing the dissipation 
rate at the truncation scale. 

Arguably, the two points above may not be significant for at- 
mospheres characterized by a very short diabatic relaxation time. 
For then the thermal damping would dominate and naturally short- 
circuit the above issues pertaining to the artificial viscosity and 
resolution. However, in some GCM simulations, the dynamically- 
relevant, intrinsic thermal relaxation time is not always short af- 
ter quasi-equilibration is reached, even above the ~1 bar level 4 
(Thrastarson (private communication)). Moreover, the instabilities 
at higher latitudes possess short growth times and are much less 
affected by short relaxation times. Additionally, there is the issue 
of transient growth, which we have not discussed in this work. The 
non-normal modes, associated with such growth, may operate on a 
much shorter time scale than the growths described in this work. 

We stress here that the high resolution runs described do not 
merely contain more fine- scale structures that presumably do not 
significantly affect the evolution. On the contrary, we have found 
that high resolution fundamentally affects the evolution. For exam- 
ple, bulk eddy heat- and momentum-fluxes differ significantly (up 
to an order of magnitude) in high- and low-resolution simulations. 
The reason for this is that vorticity anomalies (eddies) are much 
stronger, in addition to the filaments, in the high-resolution simula- 
tions; and, EKE growth is exponential in the stage of the evolution 
when these structures have emerged (Fig. 7). Hence, heat and mo- 
mentum are redistributed much more effectively. 



4 SUMMARY AND DISCUSSION 

Baroclinic instability on extrasolar planets has not been studied 
thus far. In this work we have used an advanced pseudospectral 
GCM to perform an extensive study of the stability and non-linear 
evolution of balanced jets on hot extrasolar planets. Our non-linear 
baroclinic instability calculations have been fully validated against 
previous similar calculations for the Earth (e.g. Polvani et al. 2004; 

4 It is, in general, sufficiently long below this level. 



Jablonowski & Williamson 2006). For concreteness, we have pre- 
sented here results for a model planet with physical parameters cor- 
responding to the close-in giant planet HD209458b and focused on 
the stability of high-speed (typically 1000 ms _1 ) eastward jets at 
the equator and westward jets at high latitudes. Broad jets of such 
magnitude are a common feature in current GCM simulations of 
tidally- synchronized giant planet atmospheres. 

We have derived linear growth rate and phase speed spectra, 
via standard normal mode analysis, and compared the results with 
full non-linear numerical simulations. According to our linear anal- 
ysis of the two-layer primitive equations model on the /3-plane, the 
growth rate of the instability is reduced for a jet located at low lati- 
tudes, compared with a jet located at high latitudes. Near the equa- 
tor, where the deformation length scale Ld becomes too large to 
accommodate baroclinic waves, the linear theory predicts stability. 
In general, linear analysis agrees reasonably well with the full non- 
linear calculations at the early stage of the unstable evolution, dur- 
ing the transient phase. After a long time, in simulations with high 
values of initial potential vorticity anomaly (i.e. |qf'/2f2| > 1.2, 
where q' is the anomaly), cyclones merge to form robust mono- 
lithic vortices at the poles. This is not captured by linear analysis. 

As expected, full non-linear calculations show richer behavior 
than that obtained through linear analysis. Non-linear simulations 
show that baroclinic instability occurs for all eastward jet profiles 
used in this study. In particular, broad equatorial eastward jets are 
unstable (on a time scale of ~ 20 planetary rotations), despite sta- 
bility suggested by the linear analysis. The instability takes place 
at the jet flanks, where there is still a significant vertical shear to 
satisfy the necessary condition for instability. The jet core is sta- 
ble, unlike in the jets situated at higher latitudes; this is in accor- 
dance with linear theory. Westward jets near the equator, however, 
remain stable, both at the core and the flanks. To the best of our 
knowledge this is the first time non-linear baroclinic instability has 
been studied for a broad equatorial jet in the atmosphere. In general, 
we have found westward jets to be more stable compared to their 
eastward counterparts (e.g. at midlatitude, instability timescale of 
~ 6 planetary rotations for westward jets vs. ~ 3 planetary rotations 
for eastward jets), and to require much stronger vertical shear for 
instability in the full primitive equations system. Additionally, we 
have demonstrated in this work that baroclinic instability does not 
require a solid boundary on planets, as long as there is a change 
of sign in the meridional potential vorticity gradient dqo / d<j) in the 
domain's interior. 

By performing the simulations described above with a wide 
range of horizontal resolution (from T21 to T170), we have found 
that the calculations do not converge for resolutions below T85. 
This is a somewhat stronger requirement than for Earth simula- 
tions and is primarily due to the much stronger jet amplitude on 
hot extrasolar planets (~ 1000 m s _1 , compared to ~ 50 m s _1 for 
the Earth). Furthermore, we have found that baroclinic instability 
does not occur at all if the artificial viscosity coefficient used in 
the calculation is too high. A high artificial viscosity is often used 
to stabilize numerical simulations against strong forcing in current 
studies of extrasolar planet atmospheres. Given this, baroclinic in- 
stability is unlikely to be represented in current simulations - even 
when necessary conditions for instability are satisfied. This may 
pose a serious issue in flow modelling studies of extrasolar planet 
atmospheres in which the natural diabatic relaxation time is not too 
short (i.e. greater than a few planetary rotations). 

The results presented in this paper show that baroclinic insta- 
bility is significant for understanding characteristics of hot extraso- 
lar planets which possess fast jet streams. This instability is likely 
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to play a role in weather, general circulation and large scale vari- 
ability of a few to few tens of planetary rotation periods on these 
planets (Cho 2008). For example, the process could generate large 
long-lived storms that could be observed remotely. Sharp fronts 
produced in baroclinic instability life-cycles can also act as a grav- 
ity wave source (e.g. O'Sullivan & Dunkerton 1995; Plougonven & 
Snyder 2007). Gravity waves are expected to play an important role 
in stably stratified atmospheres of hot extrasolar planets: they can 
modify the circulation through exerting accelerations (positive and 
negative) on the mean flow, as well as transporting heat vertically 
from deep regions to sensible regions and laterally from day side to 
night side (Watkins & Cho 2010). 

In this work, we only discuss adiabatic calculations. Adiabatic 
calculations are important to cleanly delineate many subtle effects 
in rotating- stratified fluid that could obscure baroclinic instability. 
They are also important as a foundation for the instability under 
forced conditions, which require careful study. The complex effects 
of forcing on the background flow itself remains to be elucidated. 
We have also focused mainly on the instability and subsequent evo- 
lution of jets in isolation and only lightly touched on the effect of 
concomitant eddies on the background flow. In summary, the full 
effect of baroclinic instability on the mean flow on hot extrasolar 
planets remains to be carefully studied. Some of the issues identi- 
fied here will be addressed in future work. 



The Charney-Green number measures the relative importance of 
the planetary vorticity gradient to the relative vorticity gradient. Us- 
ing these non-dimensionalised scales and parameters, the equations 
now read: 
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where ( r ) denotes non-dimensional variables. The vari- 
ables, G, x, 4), are the non-dimensional counterparts of 
X-> $'-) in section 2.3. 
Denoting disturbances by 

* = * exp{i/c (x -ct)} , 

where* = (^,G,x,l>) T , * = 6, X, $) T > and c G C, 
equations (Al) reduce to 

= 0, 
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APPENDIX A: NON-DIMENSIONAL STABILITY 
ANALYSIS 

Equations set (6) is made non-dimensional by introducing the 
'discretized deformation length scale', horizontal length scale, 
timescale, and height scale: 

Lb = \- \/h 2 Apao , 
Jo 

L = 7i Ld ' 



t = 



respectively. Hence, the set of equations becomes characterised 
only by the Rossby number, 



Ro = 



UoV2 



and the Charney-Green number, 



2U ' 



where 



M 



-c-j/k 2 1 

1 -c-7/P -i/(kRo) 

i/(kRo) -c-j/k 2 -i/(kRo) 

1 -ik/Ro c 



This leads to a normal mode solution fulfilling the fourth order 
characteristic equation for c : 



-4,-3/ 37 
C + C 1 



k 2 



k 6 



27 



k A 

7 



k 2 Ro 2 



&Ro 2 



Ro 2 

1_ 
k 2 



k 2 Ro 



1 



Ro 2 



1 



1 

k 2 



f 
k 4 



0. 



(A2) 



Equation (A2) is solved numerically for varying values of Ro and 
7. Fig. Al shows the results for HD209458b (cf. Fig. 2). 

Non-dimensional analysis is useful because it explicitly gives 
the dependence of stability properties on dynamically-significant 
non-dimensional numbers, such as Ro and 7. Extensive exploration 
of the solutions for a continuum of Ro and 7 values reveals that, 
as Ro and/or 7 increase, the low wavenumber cutoff for instabil- 
ity increases while the high wavenumber cutoff for instability and 
growth rate decrease slightly (Fig. Al). Fig. A2 illustrates how the 
growth rates depend on Ro when 7 is held fixed at a typical mid- 
latitude value for HD209458b. As Ro increases from Ro « 1 to 
Ro ~ 1, the growth rate decreases linearly. However, the reduction 
in growth rate is exponential as the Ro ~ 1 threshold is crossed 
and the two-layer linear analysis predicts stability for flows with 
Ro > 3.3. 

To obtain the dimensional values of the growth rate from 
Fig. Al and Fig. A2, multiply the growth rate ( k • Sm{c} ), for ex- 
ample, by UoV2/Lf> and the wavenumber (k) by V2/Lr>. The re- 
sult of Fig. Al is growth rates identical to those presented in Fig. 2. 
The dimensional phase speed is obtained similarly. 
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Figure Al. Non-dimensional growth rate [k • ^m{c}] (left) and phase speed [3fte{c}] (right) for HD209458b, as a function of non-dimensional wavenumber 
k. Curves 'HD60', 'HD45', 'HD35' and 'HD25' represent growth rates and phase speeds at = (60°, 45°, 35°, 25°) computed with (Ro= 0.76,7 = 0.14), 
(Ro= 0.76, 7 = 0.3), (Ro= 0.76, 7 = 0.51), (Ro= 0.76, 7 = 1.04). Curve 'HD45L' has been computed for HD209458b parameters at cf) = 45°, but with 
Uq = 200 m s _1 corresponding to (Ro= 0.3, 7 = 0.72). To obtain dimensional values, multiply the growth rate by U0V2/L-D and wavenumber by v^/^d- 
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Figure A2. Non-dimensional growth rate [k • Stri{c}] as a function of non- 
dimensional wavenumber k for different values of Ro. Yellow, green and 
red curves have been calculated with Ro = 10 -5 , Ro = 1 and Ro = 3.2, 
respectively. The value of 7 is held constant at 7 = 0.3. To obtain dimen- 
sional values, multiply the growth rate by UoV^/L-q and wavenumber by 
v^/Ld- 
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